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ABSTRACT 


There are several perennial issues concerning reliability estimation in Operational 
Testing and Evaluation (OT&E). Two of these include, "how should one model the 
underlying failure distribution of a continuous-time system," and, "how can a testing 
agency use information from DT in order to reduce OT resource requirements." 

In the former issue, some OT&E analysts have questioned whether or not the 
exponential failure distribution should be used in all cases for continuous-time systems, 
and have suggested the WeibuU distribution as an alternative in some instances. In the 
latter, the notion of combining DT and OT data has been an anathema to those involved in 
OT&E, however, with ever-tightening military budgets, it may be time to revisit the issue. 

First, this thesis compares the exponential and WeibuU failure distributions in terms of 
the amount of test time needed to demonstrate, to a given level of confidence, that the true 
MTTF of a sy.stem is at least as large as the minimum acceptable value, and also in terms 
of the actual confidence level associated with the lower confidence level procedure when 
the system has an increasing (or decreasing) failure rate function. Second, the thesis 
examines the behavior of an estimator for the relationship between DT and OT failure data 
using a Monte Carlo simulation. FinaUy, the thesis introduces a hierarchical Bayes 
approach for the estimation of the relationship between DT and OT failure data when a 
gamma prior distribution is assumed. 





THESIS DISCLAIMER 


The reader is cautioned that computer programs developed in this research may not 
have been exercised for all cases of interest. While effort has been made, within the time 
available, to ensure that the programs are free of computational and logic errors, they 
cannot be considered validated. Any application of these programs without additional 
verification is at the risk of the user. 
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EXECUTIVE SUMMARY 


There are several perennial issues concerning reliability estimation in Operational 
Testing and Evaluation (OT&E). One of these is, "how should one model the underlying 
failure distribution of a continuous-time system,” and another is, "how can a testing 
agency use information from Developmental Testing (DT) in order to reduce Operational 
Testing (OT) resource requirements." 

The purpose of this thesis is to (1) investigate the behavior of a statistical procedure 
used to determine the length of time needed to test a continuously operating system when 
the system under test is assumed to have an exponential time to failure distribution, but 
when in fact the system does not have an exponential time to failure distribution and rather 
may have a time to failure distribution that more closely resembles a WeibuU distribution; 
(2) to examine the behavior of a maximum likelihood estimator for the relationship 
between DT and OT failure rates using a Monte Carlo simulation; and (3) to introduce a 
hierarchical Bayes method for estimation of the relationship between the DT and OT 
failure rates. In the first part of the thesis, the parameters of the procedure are chosen 
based on the assumption of an exponential time to failure distribution. The behavior of the 
procedure is then investigated under the assumption that the time between failures are 
independent, having a WeibuU distribution with the same mean time to failure. The 
intention is to test whether or not use of the WeibuU distribution wiU result in fewer test 
hours to verify that a system under test has a mean time to failure which meets or exceeds 
a minimum acceptable value at a given level of confidence, and also to test whether or not 
the OTAs are misstating the confidence they report in their operational testing results 
when the exponential distribution is used to model a system whose true failure rate 
function is that of a WeibuU distribution. In the second part of the thesis, a maximum 



likelihood estimator is examined in terms of the consistency and variability of the OT 
failure rate derived using that estimator compared to other methods of establishing the OT 
failure rate. Finally, in the third part of the thesis, a hierarchical Bayes model is proposed 
for the ultimate purpose of being able to describe the mean and variance associated with a 
random relationship between DT and OT failure rates. 

The comparison of the exponential and Weibull failure models suggests that 
operational testing agencies can, under certain circumstances, reduce the amount of time 
needed to demonstrate that the mean time to failure (MTTF) of a continuously operating 
system is at least as large as the minimum acceptable value (MAV) if the system's true 
underlying failure distribution is Weibull with shape parameter P 1. However, the shape 
parameter is seldom, if ever, known or even estimable without expending valuable test 
resources, so it is unlikely that an operational testing agency (OTA) could actually save 
any time or resources with this methodology. Moreover, from the perspective of 
confidence levels, the analysis suggests that confidence in one's result is not sacrificed 
when the exponential distribution is used to model a system which has a non-constant 
failure rate, since it appears that use of the exponential distribution on a truly Weibull 
system seems to cause the OTA to overestimate the number of hours needed to test the 
MAV for the MTTF to a given confidence level, and also to underestimate the confidence 
it reports in its results under the incorrect assumption of an exponential distribution. The 
thesis also points out, however, that these results may not be true in general. Therefore, it 
is concluded that there is currently little to be gained by attempting to model the time 
between system failures with the Weibull distribution as previously suggested, particularly 
if very few observations of time to failure are available. 

The second analysis indicates that the maximum likelihood method for estimating a 
rate given that system's DT failure rate as well as failure data on 

vii 


system's OT failure 





previous simUar systems produces a reasonable estimate of the actual OT failure rate 
(when the OT failure rate behaves as expected). In fact, the estimate provided by the 
maximum likelihood method appears to be nearly as good as direct observation of the OT 
failure rate obtained through actual testing. Thus, if it can be shown that inferences about 
the failure rate of a system can in fact be made on the basis of previous similar systems, 
then the maximum likelihood method of predicting a new system's OT failure rate from 
previous DT failure dau has some potential utility. In particular, this OT failure rate 
prediction could be used in ways which would aid decisions as to when to begin OT as 
well as to augment follow-on OT in an effort to reduce needed testing resources. 

The third analysis is incomplete. It was hoped that the hierarchical Bayes method 
given would provide a means by which to estimate the relationship between the DT and 
OT failure rates when the relationship is assumed to be a random variable with some prior 
distribution. Regretfully, the equations derived for the maximum likelihood estimators for 
the parameters of the gamma prior have proved to be difficult to solve satisfactorily. If a 
solution to these equations were available, it could be used to predict future failure rates 
by means of a Bayesian calculation. Investigation of such a procedure is left for a later 



I. INTRODUCTION 


There are several perennial issues concerning reliability estimation in Operational Test 
and Evaluation (OT&E). In continuous systems, such as radios, radars and other items 
that operate continuously while in use, one perennial issue is "how should one model the 
underlying failure distribution of the system?" Another is "what can be done to reduce the 
amount of testing resources needed during OT&E?" 

Historically, OT&E organizations, such as Commander, Operational Test and 
Evaluation Force (COTE), have used the exponential distribution to model the behavior of 
the reliability of a continuous-time type system. The assumption of an exponential time to 
failure distribution results in a required test period of three times a given minimum 
acceptable value (MAV) for the 80 percent lower confidence limit for the mean time to 
failure (MTTF). As long as fewer than two failures occur during this time, then COTE 
will report 80 percent confidence that the true MTTF of the system is at least as large as 
the MAV (COTE, para. 204.a.(3), 1992). COTE analysts have questioned whether or not 
the exponential distribution is the best underlying failure model to use, and have 
considered investigating alternative distributions (Madson, 1993). One possible 
alternative in question is the WeibuU distribution, which contains the exponential 
distribution as a special case. The main difference between the two underlying failure 
distributions is that the exponential distribution has a constant failure rate function 
throughout the lifetime of the system, and the WeibuU distribution has a changing failure 
rate function (either increasing, constant or decreasing) during the lifetime of the system. 
The idea is that the WeibuU distribution might provide a "better fitting" model for certain 
systems which might then reduce the amount of test time needed to demonstrate, to a 
given confidence level, that the true MTTF of a system is at least as large as the MAV. 
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A third perennial issue is that of combining Developmental Testing (DT) failure data 
and Operational Testing (OT) failure data in a way that would reduce the overall amount 
of resources needed to establish the reliability of a system. For example, Gaver and 
Jacobs (1994, p. 1) describe the following problem setting: 

A test and evaluation analyst has in his possession historical (time-between-faQure 
or equivalent) data for both the...DT...and...OT...phases of a number of projects that 
are roughly comparable: they are all C3/I systems, perhaps. It is believed that there 
may be a usefully exploitable relationship between the DT and OT data. If this 
relationship could be quantified then perhaps it could be used to augment DT data 
for a new project of the same type, thereby obtaining some anticipation of the OT 
data for that new project... 

Gaver and Jacobs (1994, p. 5) go on to devise an estimator for the relationship 
between the DT data and OT data and describe its application thusly: 

The ultimate application of the above (relationship) is to (a) anticipate failure 
patterns during Operational Testing, given data on failures during Developmental 
Testing; this might aid in the decision as to when to begin OT; and (b) to 
strengthen, and reduce uncertainty of, the post-DT estimates of...(system 
reliability)...by incorporating the DT data. 

Of course, in the preceding example, it is assumed that the relationship between the 
DT and OT failure rates is constant across all "similar" systems. It is far more likely that 
the relationship follows some probability distribution. Thus, a hierarchical approach for 
the description of this relationship may be appropriate. One such approach is outlined in 
the last part of this thesis. 

The notion of combining DT and OT data has historically been an anathema to those 
involved in OT&E. The primary reason for this is that the Operational Testing Agencies 
(OTAs) have, rightly so, desired to be completely independent of the Developmental 
Agencies (DAs). However, with the tight budgets facing all military organizations of late, 
it may be time to re-visit the idea of using previous DT data to supplement OT data in an 
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effort to reduce the amount of testing resources needed to establish a system's reliability in 
OT&E. 

The fust goal of this thesis is to compare, through a sensitivity analysis, the 
exponential and Weibull failure distributions in terms of the amount of test time needed to 
demonstrate, to a given level of confidence, that the true MTTF of a system is at least as 
large as the MAV, and to find the actual confidence level associated with the lower 
confidence limit procedure when the system has an increasing (or decreasing) failure rate 
function. The former tests whether or not use of the Weibull distribution wUl result in 
fewer test hours to establish reliability with a given confidence level. The latter tests 
whether or not the OTA is overstating (or understating) the confidence it reports in its 
results when the exponential distribution is assumed but the system's true failure rate 
follows a WeibuO distribution. 

The second goal of this thesis is to examine the behavior of the estimator for the 
relationship between DT and OT data devised by Gaver and Jacobs (1994) using a Monte 
Carlo simulation. The behavior to be examined includes the consistency and variability of 
the OT failure rate derived using this estimator compared to those of other methods of 
establishing the OT failure rate, including direct observation of the OT data and a 
weighted average between direct observation and the Gaver-Jacobs estimate of the OT 
failure rate. 

The third goal of this thesis is to introduce a hierarchical Bayesian approach for the 
estimation of the relationship between DT and OT data using a gamma prior distribution. 

Note that the methods used for borrowing DT information to assist in predicting OT 
information, or actual later field experience, can also be used in other situations. For 
example, suppose an upgrade of an existing system is to be evaluated. Then actual data 
from the existing system may be scaled to the new system in a manner analogous to that 
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already described for DT-to-OT scaling. Note also that the general procedure described 
can be adapted to other operational parameters. An example might be an effectiveness 
parameter such as projectile dispersion around an aim point, technically speaking, a 
variance parameter. Under current conditions of declining resources such an approach 
could be of considerable value in enhancing the accuracy of many types of test data. 
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n. COMPARISON OF EXPONENTIAL AND WEIBULL DISTRIBUTIONS 


A. EXPONENTIAL MODEL 

One of the simplest and most popular ways to model the distribution of the random 
time to failure of a continuous time system is with the exponential distribution. Among 
the reasons for the popularity of the exponential distribution are its simplicity, ease of use, 
and the memoryless property which says that the probability of failure within the next t 
time units is independent of the frequency of any previous failures. 

Let T = random time to failure of a system, 
t = length of time system is in operation, 

7 = failure rate. 

The time to failure, T, of a system follows an exponential failure distribution if the 
probability density function (PDF) is 


f(t) = 7e-"‘,t>0 

and if its cumulative distribution function (CDF) is 


F(t) = l-e-^',t>0 

= 0,t<0 


If the times between system failures are independent and have the exponential 
distribution specified above, then the system's mean time to failure (MTTF) is 


MTTF = E[T] = - 


and its failure rate function is 


(2.4) 





This distribution is used by the analysts and Operational Test Directors (OTDs) at 
COTF primarily to determine the length of time needed to test a given continuous time 
system in order to demonstrate that the true M l IT of the system is greater than or equal 
to the MAV for the MTTF with approximately known confidence. The basic assumption 
associated with the exponential distribution is that the system under scrutiny exhibits 
steady state performance; that is, the time between failures are independent and have a 
common exponential distribution. In other words, it is assumed that the failure rate 
function, h(t), is constant over time: there is no "reliability growth," or degradation. 

COTF's usual policy is to conduct sequential tests. That is, they will test one system 
for a pre-determined length of time. If it fails before the time on test expires then they will 
repair the system and continue the test. If the repaired system fails before the remaining 
time on test then the system fails the test, else it passes. In other words, if there are zero 
or one failures during the time on test, then the system passes (COTF, para. 602.a, 1992). 

Let a = confidence level, 

e = MAV for the MTTF, 

T2 = time to 2nd failure, 

N(t) = number of failures to occur in [0, t]. 

The time, t, needed for 100(a) percent confidence, with no more than one failure 
allowed during testing, required to demonstrate that the system meets the MAV for the 
MTTF, 0, is determined as follows: 


l-a>P(Tj>t) = P(N(t)<l) 
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(2.5) 




Standard procedure at COTF is to set the confidence level at 80 percent and to allow 
for a maximum of one failure to occur during the operational test. Solution of equation 
(2.5) with a = 0.80 gives a time to test, t = 36.i In other words, COTF will plan to test a 
continuous time system for a period of time which is three times greater than the MAV for 
the MTTF, and if the system has fewer than two failures in this time, then the system 
passes the test at the 80 percent confidence level. Note that other combinations of test 
times and acceptable numbers of failures can yield the same confidence, but the above has 
been chosen as a practical compromise. 

For example, if the MAV for the MTTF is 1000 hours, COTF will test the system for 
3000 hours and allow at most one failure to occur. If the system passes this test, COTF 
will report that it is 80 percent confident that the system's true MTTF meets or exceeds 
the MAV. 

Implicit in the preceding model, however, is the notion that if the expected value of 
the time between failures for one system is greater than or equal to the expected value of 
the time between failures for a second system, then the probability of having J or fewer 
failures in a renewal process observed for a fixed time for the first system is at least as 
large as that for the second system given the same fixed time. 

Let X = time to 1st failure of system 1, X - exp(7,), 

Y = time to 1st failure of system 2, Y - exp(72). 7i ^ 72- 

Proposition 1: If E[X] > E[Y] then P(Nx(t) < 1) > P(Ny(t) < 1) V t. 

Proof; 7, < 72 => W = l/7i ^ E[Y] = 1/y^ 

7i < 72 => Fx(t) = e-’'‘' > FY(t) = e“^‘' Vt>0 

^ An interesting, and simpler, alternative method for determining t exists, involving the chi-square 
distribution. See Appendix A for details of the alteiimlive method. 
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Let 


Nx(t) = the number of renewals that occur in a renewal 
process having inter-arrival distribution F^lt), 

Nylt) = the number of renewals that occur in a renewal 
process having inter-arrival distribution Fylt). 

then Fx (t) > Fy (t) => Nyd) V t > 0 (Ross, p.257, 1983), 

which means that P(Nx(t) < 1) > PCN^Ct) < 1). Q.E£>. 

Note that Proposition 1 may not be true for aU eligible distributions. However, the 
fact that it is true for the exponential distribution enables one to use the COTF procedure. 

B. ON CONFIDENCE LEVELS 

One might ask at this point how COTF defines "level of confidence." As used in 
equation (2.5), the confidence level is actually the "producer's risk." When COTF desires 
to test a system’s MTTF with 80 percent confidence, this means they are willing to assume 
a 20 percent (or smaller) risk of incorrectly recommending a system for production based 
on the results of their test. 

Let Nx(t) and NY(t) denote, respectively, the number of renewals to occur in two 
renewal processes having distribution Fx(t) and FY(t). Ross (p. 257, 1983) shows that if 
Fx(t) > FY(t) then P(Nx(t) < 1) > PINyIO < 1). Since Fx(t) > FY(t) inplies that E[X] > 
E[Y], the result shows that under this additional assumption, if the true MTTF of the 
system is less than the MAV, 0, and the desired level of confidence is a, then Ross's result 
implies that P(N(t) < 1) < 1 - a and P(N(t) > 1) > a. 

In other words, if COTF conducted the same test many times on a system with 
exponential times between failures, and if the true MTTF were actually less than the 
MAV, 0, then the observed time for the second failure to occur, T2, would fall below 30 
more than 80 percent of the time (in which case COTF would correctly reject the system). 
Consequently, less than 20 percent of the time the observed time for the second failure to 



occur, T2, would lie above 36 (in which case they would incorrectly accept the system). 
Symbolically, COTF desires; 

P(accept the system I true MTTF < MAV) < .20. 

There is an interesting implication to this requirement. In order to have any 
appreciable chance of passing a test set up according to COTF standards, the DA must 
build its system such that its true MTTF is much greater than the MAV (KeUer, 1993). If 
the system's MTTF is equal to the MAV, then it has only a 20 percent chance of pa.ssing 
the test. 

Note carefully: the above does nsu say that if the test leads to acceptance, i.e., that 
fewer than two failures are observed in three times the MAV, then the probability is 80 
percent that the true MTTF is greater than the MAV. This statement is specifically not 
allowed in the present framework. 

C. WEIBULL MODEL 

An alternative, and more general, way to model the probability of failure of a 
continuous system is with the WeibuU distribution. The Weibull distribution is useful 
when the system under scmtiny displays either an increasing or decreasing failure rate over 
time. An example of a system with an increasing failure rate function would be one that 
"wears out" over time. A system that exhibits "wear out" would tend to have a higher 
failure rate in old age than an identical system which is younger. An example of a system 
with a decreasing failure rate function would be one that "wears in" over time. Certain 
electronic devices may exhibit "wear in" tendencies, at least initially. 
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Let T = random time to failure of a system, 
t = length of time system is in operation, 

X = scale parameter (characteristic life), 

P = shape parameter. 

The time to failure, T, of a system follows a WeibuU failure distribution if the 
probability density function (PDF) is 


f (t) = p)t^t^-‘e-<^'", t > 0, P > 0, >. > 0 


( 2 . 6 ) 


and if its cumulative distribution function (CDF) is 


F(t) = l-e-'^>’,t>0,P>0,X>0 
= 0,t<0. 


(2.7) 


If the times between system failures are independent and have the WeibuU distribution 
specified above, then the mean time to failure (MTTF) of the system is 


d its failure rate function is 




(2.8) 


(2.9) 


This distribution has been considered as a replacement for the exponential distribution 
with the hope that it might more precisely model many of the systems under test by COTF. 
In this model, X is analogous to the exponential failure rate, y. The shape parameter, P, 
determines whether the system's failure rate is increasing, decreasing, or constant over 
time. When P is less than one, the failure rate function is decreasing, when P is greater 
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than one, the failure rate function is increasing, and when P is equal to one, the failure rate 
function is constant. Note that when p is equal to one, the WeibuU distribution is the 
exponential distribution. 

There is a useful relationship between the exponential MTTF, 0 = I/7, and the 
WeibuU parameters ^ and P (when the exponential and WeibuU MTTFs are equal): 



This relationship wiU be useful when comparing the performance of the two distributions. 

Let a = confidence level, 

T2 = time to 2nd faUure, 

N(t) = number of failures to occur in [0, t]. 

The time, t, needed for 100(a) percent confidence, with no more than one failure 
allowed during testing, required to demonstrate that the system meets the MAV for the 
MTTF, 6, may be determined, using an adaptation of the test for the exponential 
distribution, equation (2.5), to the WeibuU distribution as foUows: 

1 - a > P(T2 > t) = P(N(t) < 1) = P(N(t) = 0) + P(N(t) = 1) 

P(N(t) = 0) = 

P(N(t) = 1) = JFt, (t-x)fT.^ (x)dx = PX^ j x^‘e'‘'^-^>“-^‘’^>“'dx 
1 - a > e-'"**' + P>.^ j x^-‘e-'‘’‘-’">'‘^''‘’"''dx (2.11) 
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Solution of equation (2.11), with a, P, and X specified gives the appropriate time to 
test, t. Note, however, that there is no closed form solution for t in this inequality. It 
must be solved by numerically integrating for a number of t-values. 

Implicit in the preceding model, again, is the notion that if the expected value of the 
time between failures for one system is greater than or equal to the expected value of the 
time between failures for a second system, then the probability of having J or fewer 
failures in a renewal process observed for a fixed time for the first system is at least as 
large as that for the second system given the same fixed time. 

Let X = time to 1st failure of system 1, X - Weibull( p, X,), 

Y = time to 1st failure of system 2, Y - WeibuUf p, X2), X, < , P > 0. 


Proposition 2: If E[X] > E[Y] then P(Nx(t) < 1) > P(NY(t) < 1) V t. 


Proof; X, < Xj 


n.M.Ervi.lM 


X, <X2 => Fx(t) = e"<’-‘'" > FY(t) = e-‘*-"‘’' Vt>0,P>0 

Let Nx(t) = the number of renewals that occur in a renewal 
process having inter arrival distribution Fx(t), 

NyCt) = the number of renewals that occur in a renewal 
process having inter arrival distribution Fyft). 

then Fx(t)> FY(t) => Nx(t) ^NyCt) V t > 0 (Ross, p.257, 1983) 

which means that P(Nx(t) < 1) > P(NY(t) < 1). Q.EJ). 

D. TEST TIME NEEDED TO DEMONSTRATE A THRESHOLD 


Suppose COTF is testing a new radar system whose MAV for MTTF is 1000 hours. 
Suppose further, for the time being, that the radar system's true MTTF is in fact 1000 
hours, i.e., y = 0.001. Using the assumption of a constant failure rate function, and the 
thumb rule of 80 percent confidence with no more than one failure, COTF would then test 


12 



the system for 3000 hours and allow at most one failure to occur. If the system passes this 
test, COTF will report that, with 80 percent confidence, the radar system MTTF meets or 
exceeds the MAV. 

Now suppose COTF wishes to model this system’s behavior with a WeibuU 
distribution. Several choices for the shape parameter, p, are available. Assume that, 
depending on the platform on which the radar system will be deployed, the shape 
parameter may vary anywhere between 0.5 and 3.0. The relationship derived in equation 
(2.10) can now be used to match the MTTF of the exponential distribution, I/7, to the 
MTTF of a Weibull distribution given various shape parameters, p, as seen in Table 1. 


TABLE 1 

RELATIONSHIP OF EXPONENTIAL FAILURE RATE, 7, TO WEIBULL SCALE 
PARAMETER, X, TO MATCH MTTF's FOR DIFFERENT WEIBULL SHAPE 


PARAMETERS, P 


Exponential 
failure rate, 7 

Weibull 

shape parameter, p 

ni/p) 

Weibull 

scale parameter, X 

0.001 

0.5 

1 

0.002 

0.001 

1 

1 

0.001 

0.001 

2 

1.7725 

0.000886 

0.001 

3 

2.6789 

0.000893 


Now it is possible to use equation (2.11) to derive the appropriate test times for the 
preceding values of P as shown in Table 2, again assuming a confidence level of 80 
percent, and allowing for at most one failure to occur.2 


^The author solved equation (2.11) using MathCad version 3.1, varying t in one hour increments until the 
inequahty was satisfied. 
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TABLE 2 

TEST TIMES ASSOCIATED WITH VARIOUS WEIBULL 
_ SHAPE PARAMETERS, p _ 


p 

X 

t (hrs) 

0.5 

0.002 

2950 

1 

0.001 

3000 

2 

0.000886 

2613 

3 

0.000893 

2435 


This analysis suggests that, everything else being equal, the assumption of an 
exponential failure distribution (|3 = 1) generates the longest test time to demonstrate that 
the system meets or exceeds the MAV for MTTF at an 80 percent confidence level. In 
other words, if it is reasonable to assume that the system under test follows a WeibuU time 
to failure distribution, and if the OTA can identify a reasonable P value for the system 
(where P 1), then by modeling the system with a WeibuU failure distribution the OTA 
may be able to reduce the required number of hours needed to demonstrate that the 
system meets or exceeds the MAV for MTTF for a WeibuU distribution with the same 
shape parameter while maintaining the same level of confidence in its result 

E. USE OF EXPONENTIAL MODEL ON A WEIBULL PROBLEM 

Now suppose COTF uses its traditional rule for determining time needed to test the 
radar system. Earlier it was shown that for an exponential failure distribution with a MAV 
for MTTF of 1000 hours, it was necessary to test the radar system for 3000 hours, 
aUowing for at most one failure, to obtain 80 percent confidence that the true MTTF met 
or exceeded the MAV. Suppose it is later discovered that the true underlying failure 
distribution was actually WeibuU with shape parameter P 1. What can then be said 
about the actual level of confidence achieved from the test? Is it still 80 percent? 

One way to explore this is to hold the test time constant at 3000 hours, and determine 
the probability of passing the test ( 1 - a ) and also the level of confidence ( a ) achieved 
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when the underlying failure distribution is WeibuU. This can be accomplished using 
equation (2.11) by fixing the time to test, t, at 3000 hours, specifying |3 and X, and then 
solving for a. Results are shown in Table 3.3 

This analysis seems to suggest that no confidence is lost by assuming an exponential 
failure distribution when the true underlying failure distribution is WeibuU with parameter 
P > 1 . This makes intuitive sense. As was suggested back in Table 2, for any P > 1, the 
required test time is reduced, thus if the system is tested longer than necessary to achieve 
80 percent confidence, as occurs in Table 3 when P > 1, then it naturally follows that the 
actual confidence level in the result of the test should increase. 


TABLE 3 

ACTUAL CONHDENCE LEVELS ACHIEVED AS A FUNCTION OF 
THE SHAPE PARAMETER, p, WHEN EXPONENTIAL 


p 

X 

t(hrs) 

P(pass test) 
(1-a) 

Conf. level 
(a) 

0.5 

0.002 

3000 

0.196 

0.804 

1 

0.001 

3000 

0.200 

0.800 

2 

0.000886 

3000 

0.097 

0.903 

3 

0.000893 

3000 

0.028 

0.972 


Unfortunately, the results of this analysis may not hold under aU circumstances. The 
reason for this is because Propositions 1 and 2 no longer apply. In other words, the fact 
that E[X] > E[Y] no longer implies that P(Nx(t) < 1) > PCNyCt) < 1). 

Let X = time to 1st failure of system 1, X ~ exponentially), 

Y = time to 1st failure of system 2, Y - Weibull( P, X). 


^Results obtained by solving equation (2.11) using MathCad 3.1, varying a until the inequality was met. 
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Assume; E[X] = - > E[Y] = 

7 P>- 

Unfortunately — < —does not imply that 
7 pA 

Fx(t) = e'^^'' > FY(t) = e’<*-‘’' Vt>0,P>0. 

Therefore, P(Nx(t)<l) > P(NY(t)<l) cannot be shown using the proof given for 
Propositions 1 and 2. 

Thus, the reader is cautioned not to generalize the result obtained in the previous 
analysis. It is unclear how the renewal process with exponentially distributed inter arrival 
times compares to the renewal process with WeibuU distributed inter arrival times when 
their mean arrival times are equated. 

F. CONCLUSIONS 

The preceding analyses suggest that the OTA can in fact reduce the amount of time 
needed to test a system for mean time to failure if the system under test follows a WeibuU 
faUure distribution with shape parameter p > 1, and if this shape parameter is known or 
can be estimated. Unfortunately, the shape parameter p is seldom, if ever, known, and it is 
also difficult to estimate without prior testing, which tends to make moot the notion of 
using the WeibuU faUure distribution to reduce needed test resources.^ So, fi'om this 
perspective, one can conclude that there is probably Uttle utility in the OTA attempting to 
model system reUabUity using the WeibuU distribution as previously suggested. 

Moreover, from the perspective of confidence levels, it remains unclear how the use 
of the test based on the exponential distribution, when applied to a WeibuU problem, 

'‘The Gaver-Jacobs method for "borrowing information" from previous similar systems may be useful in 
this regard. Details of the Gaver-Jacobs method is given in the next section. 
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affects the actual confidence levels achieved. Analysis suggests that, for purposes of 
making inferences about the mean time to failure, it might be perfectly safe for the OTA to 
use the exponential failure distribution to model a system that would perhaps more 
appropriately be modeled using the Weibull distribution, since, as suggested in Tables 2 
and 3, modeling a system which has a Weibull failure rate function with the exponential 
distribution seems to cause the OTA to overestimate the number of hours to test, which in 
turn seems to provide more confidence in the result of the test. It will be desirable to 
investigate conditions under which this conservatism property is generally true. 

Finally, it should be remembered that OTAs such as COTF tend to perform testing at 
the system level in which there are numerous competing failure modes. It is well- 
established in the literature (e.g., NAVORDSYSCOM, p. 3-25, 1971) that as the number 
of competing failure modes increases, the occurrence of any one failure tends to occur at a 
constant rate, regardless of the individual underlying failure distributions. So, despite the 
warning concerning modeling a Weibull system with an exponential time between failure 
distribution, for purposes of making inferences about the mean time to failure, with larger 
systems it is not necessarily inappropriate to assume an exponential failure distribution 
with a constant failure rate function. 
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m. ON COMBINING DT AND OT FAILURE DATA 


A. BACKGROUND AND PROBLEM DESCRIPTION 

It has recently been postulated that there may be a useful relationship between 
Developmental Testing (DT) and Operational Testing (OT) time-between-failure data for 
many systems. In fact it has been noted anecdotally by some researchers that the OT 
failure rate tends to be roughly four times the DT failure rate for a given system.^ Gaver 
and Jacobs (1994) note that if this relationship could be quantified, then it could be used 
to anticipate the OT failure rate of a new, similar system based on the DT failure rate from 
that system. 

1. Example Problem Setting 

An OT&E analyst has in his possession historical time-between-failure data for 
both the DT and OT phases of ten projects that are similar, e.g., radar systems. The 
analyst believes a relationship exists between the DT and OT failure rates such that 

A.„(i) = ia,(i) (3.1) 

where X„(i) represents the OT failure rate of system i, X,j(i) represents the DT 
failure rate of system i, and K represents a constant of proportionality^ and where the 
times between failures are independent, having an exponential distribution with means 
l/^„(i) and 1/A.^ (i), respectively. The analyst further believes that if he can determine the 

^Interview between D. P. Gaver, Professor of Operations Research, Naval Postgraduate School, and the 
author. March 1994. 

^The author recognizes that this "constant" of proportionality is actually a random variable with a 
distribution of its own. Refer to section E of this chapter for a Bayesian method of estimation of K. 
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value of K, then he can use this value on an 11th similar system to anticipate the OT 
failure rate of System 11 given System 11's DT failure data. 

B. THE GAVER-JACOBS MODEL 

Suppose T(](i;j) is a random variable representing the time between the (j-l)st and the 
jth time to failure in DT of system i where failures occur according to a Poisson Process 
such that 


P{T,(i,j)>tJ = e-^-'‘>‘^ , t, >0. (3.2) 

Suppose further that a similar relationship holds for To(i;j), the random variable 
representing the time between the (j-l)st and the jth time to failure in OT of system i 

PK(i.j)>t„} = e-^-''>'- , t„>0. (3.3) 

Then N£j(t£j(i)) is the number of failures to occur during DT for a test of length t£j(i) 
fixed in advance, and has a Poisson distribution, i.e.: 

P{N,(t,(i)) = n} = ^ n = 0,l,.... (3.4) 

n! 

An equivalent distribution holds for OT. 

Assume that for each system i, there are n^jli) DT failures in a test of length tcj(i), and 
no(i) OT failures in a test of length to(i). Gaver and Jacobs (p. 3, 1994) have devised a 
maximum likelihood estimator for the constant of proportionality K as follows: 

The likelihood function is 


L(L,L;data) = ne"^“‘ 




n,(i)! 


n„(i)! 


(3.5) 
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where I denotes the total number of similar systems. Now assume = K?Lj(i); then, 
up to irrelevant constants, the likelihood of K and the development parameters is 


L(K,A.j(i);data) = e ‘ ' 


‘.(i) I 

nA.,(i)”''‘>K''”%(i)"« 


(3.6) 


where n„ = ^ n„ (i), the total number of OT failures across all systems. Focusing on the 
parameter K above, take the natural logarithms of both sides of equation (3.6) to get 


^(K, A.,(i);data) = n> K -X^d(i)(td(i) + Kt„(i)) 
then differentiate with respect to K to get 


(3.7) 


then finally 
estimator 


aK 

equation (3.8) equal to zero 


(3-8) 

and solve for K to get the maximum likelihood 


K = ^ (3.9) 

A.dt„ 

where A-^t^ = ^ (i)t„(i), and to(i) is the length of time over which system i is subjected 

to test in OT. Unfortunately A^ (i) is unknown, but may be estimated as follows: 

A,(i)=i^i^. (3.10) 

td(i) 

The properties of this estimator for K are as yet unknown and will be investigated 
shortly. Meanwhile, notice that, assuming K is well-behaved, the analyst now has a 
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mechanism for relating past DT failure data on a system to future OT failure data on the 
same system given historical DT and OT failure data from similar systems. Specifically 

X„(i) = id,(i). (3.11) 

C. ANALYSIS OF THE GAVER-JACOBS MODEL 

The test of the performance of the estimator, K, and the follow-on estimate ^„(i), 
was conducted using a Monte Carlo simulation. The idea was to simulate ten "similar" 
systems (Systems 1 through 10) whose actual OT failure rates were four times their actual 
DT failure rates. Each system was tested for a length of time, t, such that the mean 
number of failures in (0, t], A.t, was equal to 4.0 for each system.^ Associated with each 
system was a random number seed which was used with a Poisson random number 
generator, written by the author in Turbo Pascal 6.0, the code for which is included in 
Appendix B. A summary of the initial data is given in Table 4. 

1. Estimation of K. 

Next, 1000 separate observations of the number of OT failures and DT failures 
observed over the appropriate test time for each system were generated using the random 
seed indicated in Table 4. Then, K was derived using equations (3.9) and (3.10) for each 
of the 1000 repetitions yielding 1000 separate observations of K. Appropriate statistics 
were gathered on K. Table 5 summarizes these statistics and Figure 1 shows a histogram 
of the observations of K. 


^This number was chosen arbitrarily, but was kept constant to simulate consistent choice of test duration. 
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TABLE 4 

INITIAL DATA ON TEN "SIMILAR" SYSTEMS IN WHICH 
ACTUAL K = 4.0 AND ACTUAL POISSON FAILURE RATE A,t = 4.0 


System 

number 

Actual DT 
failure 

rate 

DT test 
time 
(hours) 

DT 

Random 

seed 

Actual OT 
failure 

rate 

OT test 
time 
(hours) 

OT 

Random 

seed 

1 

0.0002 

20,000 

10 

0.0008 

5,000 

110 

2 

0.0004 

10,000 

20 

0.0016 

2,500 

120 

3 

0.0006 

6,666.67 

30 

0.0024 

1,666.67 

130 

4 

0.0008 

5,000 

40 

0.0032 

1,250 

140 

5 

0.001 

4,000 

50 

0.004 

1,000 

150 

6 

0.002 

2,000 

60 

0008 

500 

160 

7 

0.004 

1,000 

70 

0.016 

250 

170 

8 

0.006 ^ 

666.67 

80 

0.024 

166.67 

180 

9 

0.008 

500 

90 

0.032 

125 

190 

10 

0.01 

400 

100 

0.04 

100 

200 


TABLES 

SUMMARY STATISTICS ON 1000 
OBSERVATIONS OF K 


Mean 

4.040 

Std Error of the Mean 

0.030 

Median 

4.000 

Std Dev of the simulated values of K 

0.939 

Variance of the simulated values of K 

0.882 

Kurtosis 

0.677 

Skewness 

0.675 

Range 

5.779 

Minimum 

1.913 

Maximum 

7.692 

Sum 

4039.955 

Count 

1000 
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Figure 1 - Histogram of 1000 observations of estimate of K. 

2. Estimation of the OT failure rate. 

While information about the sampling distribution of K is interesting, the real 
focus of this study is on how well the analyst can predict the OT failure rate. Three 
methods for estimating the OT failure rate are being considered in this study: 

1. Prediction using the relationship (i) = (i); 

2. Direct observation of the OT failure rate from OT data; and 

3. Weighted average of the two previous methods. 

In each of the three methods, the following are calculated: 

1. The estimate for the OT failure rate; 

2. The variance of the estimate, derived appropriately from each method; 

3. The standard error (square root of the variance) of the estimate; and 

4. The mean-squared error, which is the sum of the square of the average of the biases 
and the average of the variances over all replications. 
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Figure 2, below, describes one replication of the OT failure rate estimation 


simulation model. 



Figure 2 - One replication of the OT failure rate estimation simulation model. 


a. Prediction using the estimator for K. 

The analyst’s prediction of the future OT failure rate for an 11th system, 
using the estimator for K, would be 


X„(ll) = ia,(ll). (3.12) 

The associated variance would be computed as follows: 

Var[X„(ll)] =Var[id.,(ll)] 

= E[(^,(11)-E[K]E[X,(11)])2] 

= E[K2^, (11)2 - (11)E[K]E[X, (11)] + E[K]2e[^, (11)]2] 

= E[K2]E[^, (11)2] - E[K]2e[X, (11)]2] 

= (Var[K] + E[K]2 )(Var[^, (11)]+ E[i, (11)]2 ) - E[K]2e[^, (11)]2. 
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simplifying gives; 

Var[^Jl I)] = Var[K]Var[i,(ll)] + Var[K]E[X,(ll)]2 ) + Var[X.,(ll)]E[K]2 . (3.13) 

Expressions are needed for Var[K] and 'Var[Xj(ll)]. An approximation to 
Var[K] can be derived through maximum likelihood estimation using observed Fisher 
information as follows (Cox and Hinkley, p.297, 1974): 

Var[K] = —-. 

Taking the second derivative with respect to K of equation (3.7) gives 
dV n„ 

die" K"' 


(3.14) 


(3.15) 






Var[^j(ll)] can be derived as follows by noting that ^^(i) is simply a Poisson random 




(t.(i)) (t,(i)) 
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The standard error of (11) is then simply the square root of Var[^j (11)]. 

An example prediction can now be performed. Suppose, using data from the 
fust simulation run, equation (3.9) gives K = 4.293 and equation (3.10) gives A.j(ll) = 
0.0003. Then 

X„(ll) = ia,(ll) = (4.293)(0.0003) = 0.00129. 


Note that per the original assumption, K = 4.0, so we would expect a system whose DT 
failure rate is 0.0003 to have an OT failure rate of 0.0012. So, at first glance this appears 
to be a relatively good estimate of the true situation. 

The variance and standard error of X„(ll) is solved using equation (3.13) as 

follows: 


E[K] = K =4.293 
E[X,(ll)] = Xj(ll) =0.0003 
Var[K] = 




00.25)' 


Var[)t011)] = 


(t,(ll))' (20,000)' 


finally 

Var[^„(ll)] = Var[ia,(ll)] = (0.4188)(1.5xl0-8) + (1.5xl0-8)(4.293)2 
+ (0.4188)(0.0003)2 = 32 x lO'S. 


And the standard error of ^„(11) is V32xl0'* = 5.66 x lO"*. Table 6 summarizes the 
results of this prediction. 
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TABLE 6 

SUMMARY OF RESULTS OF PREDICTION OF OT FAILURE RATE 
ON SYSTEM 11 USING ESTIMATE OF K BASED ON FIRST REPLICATION 


Actual Failure Rate (11) 

0.0012 

Estimated Failure Rate (11) 

0.00129 

Variance[^„(ll)] 

32 X 10-8 

Standard EiTor[i„ (11)] 

5.66 X 10-4 

(11) ± 1 Standard Error 

(0.0007,0.0019) 

^„(11) ±2 Standard Errors 

(0.0002, 0.0024) 


b. Direct observation from OT failure data. 

Of course, the analyst could, if he desired, come up with a value for 
by observing the empirical OT failure rate as usual. In this situation, the analyst would 
note the number of failures which occur, no(ll), during the time on test, to(ll), and 
derive the OT failure rate as follows: 




n„(ll) 

todl)’ 


(3.18) 


For example, suppose that System 11 is operationally tested for a period of 
5,000 hours and experiences 6 failures during that time. Then the OT failure rate 
observed through operational testing would be 


Using the same reasoning as in equation (3.17), the variance of the observed 
value is calculated as 
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Var(N„(ll)) 


Var[X:(ll)] = Var 


No(ll) ] = f 1 T 


=-^ = 24x10'*. 

(5000)" 

And the standard error of X'(ll) is ■\/24xl0“* = 4.9 x lO""^. Table 7 summarizes the 
results of these calculations. 


TABLE? 

SUMMARY OF RESULTS OF DIRECT OBSERVATION OF OT FAILURE RATE 
_ON SYSTEM 11 FROM OT FAILURE DATA_ 


Actual Failure Rate (11) 

0.0012 

Estimated Failure Rate X' (11) 

0.0012 

Variance[X'(ll)] 

24 X 10-8 

Standard Error[X'(11)] 

4.9 X10-4 

X'(ll) ± 1 Standard Error 

(0.0007,0.0017) 

X' (11) ± 2 Standard Errors 

(0.0002,0.0022) 


Unfortunately, this method would not save the analyst any time since he 
would have to run the entire operational test procedure to derive a meaningful number. 
Moreover, he would not be taking advantage of the available DT failure data. 
c. Weighted average of predicted value and direct observation. 

Suppose now that the analyst runs the entire operational test procedure as 
above, but also desires to use the available DT failure data in an effort to improve his 
estimate for the OT failure rate. Assuming the analyst has the data from the 10 other 
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previous systems with which to derive the prediction, A,^(ll) = lOij(ll), he can then 
calculate a weighted average of the two estimates. 

An estimator for the mean of a linear combination of two random variables, 
X and Y, with common mean m and different variances, and Oy can be constructed as 
follows: 

Let m=aX + bY, 0<a,b^l, a + b=L 

Then E[m] = am + bm = m(a + b) = m(l) = m, 

^(a) =E[(m-m)2] = E[(aX + (l-a)Y-m)2] 

= E[(a(X-m) + (l-a)(Y-m))2] 

= E[a2(X-m)2 + 2a(l-a)(X-m)(Y-m) + (l-a)2(Y-m)2] 

= a2Var[X] + (1 - a)2Var[Y]. 

Differentiating with respect to a and setting equal to zero gives 

= 2aVar[X] - 2(1 - a)Var[Y] = 0 
9a 

which yields 


a[Var[X] + Var[Y]] = Var[Y] 


or 


Var[Y] , 
Var[X] + Var[Y] + 


Similarly 

b=_MX]_^ /of 

Var[X] + Var[Y] + 
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Substitution then yields 


Therefore 


_ Xx %y 


estimator for the OT failure 


rate, >.„(i),^ 


(3.19) 


be constructed as follows: 


iad(i) I ^;(i) 

Var[id,(i)] Var[X;(i)] 
1 ^ 1 ' 

Var[id,(i)] Var[l'(i)] 


where KXj (i) represents the predicted value of the OT failure rate and A.'(i) denotes the 
empirical OT failure rate observed during operational testing. 

For example, the values obtained previously for I^j(ll), the predicted 
value, and X,'(ll), the observed value, were 0.00129 and 0.0012 respectively. It was 
shown earlier that the estimated variance for the predicted value was Var[i^j(ll)] = 32 
X 10‘8, and the variance of the observed value was Var[X'(ll)] = 24 x lO"^. Inserting 
these numbers in equation (3.20) gives 

12.9x10^ ^ 12xl0~^ 

X„(ll)= 32 xLQ-* - 24x10'^ =12.39x10^. 

32x10'* "^24x10-* 

The variance associated with this estimate can be constmcted as follows: 

Let Var[m] = Var[aX + bY] = a2Var[X] + b2var[Y]. 
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Substitution yields 


(/gJPx+UJ'gy 1 
(Xi+Xj Xx'^/^Y 


Therefore, the estimated variance of A.^ (11) is computed as follows: 
Var[X„ (11)] =-j--j-=-j-i-j- 


Var[ia,(ll)] Var[l;(ll)] 32x10-* 24x10^ 

And the standard error of ^„(11) is VlS.TxlO"* = 3.7 x 10"^. Table 8 summarizes the 
results of these calculations. 


TABLE 8 

SUMMARY OF RESULTS OF WEIGHTED AVERAGE OF 
PREDICTED VALUE AND DIRECT OBSERVATION OF OT FAILURE RATE 


Actual Failure Rate (11) 

0.0012 

Estimated Failure Rate (11) 

0.00124 

Variance[X„(ll)] 

13.7 X 10-8 

Standard Error[X„ (11)] 

3.7 X 10-4 

(11) ± 1 Standard Error 

(0.0009, 0.0016) 

X„ (11) ±2 Standard Errors 

(0.0005, 0.0020) 


3. Comparison of the three methods for estimating the OT failure rate. 

Table 9 summarizes the results obtained for each of the three methods of deriving 
the OT failure rate for one replication of the simulation. They are arguably quite similar. 
So to answer the question of which method is best, further analysis is required. 
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TABLE 9 

COMPARISON OF THE THREE METHODS FOR OBTAINING 
_THE OT FAILURE RATE_ 


Statistic 

Prediction Method 

Direct Observation 
Method 

Weighted Average 
Method 

Actual Value 

0.0012 

0.0012 

0.0012 

Estimated Value 

0.00129 

0.0012 

0.00124 

Variance 

32 X 10-8 

24 X 10-8 

13.7 X 10-8 

Standard Error 

5.66 X 10-4 

4.9 X 10-4 

3.7 X 10-4 

Estimate ± 1 S.E. 

(0.0007,0.0019) 

(0.0007, 0.0017 ) 

(0.0009, 0.0016) 

Estimate ± 2 S.E.s 

(0.0002, 0.0024) 

(0.0002,0.0022) 

(0.0005, 0.0020) 


Since the results of the three methods discussed previously were so similar, a 
look at the asymptotic behavior of the three methods is in order. To facilitate this 
analysis, four new "similar" systems were created, systems 11, 12, 13 and 14, each with 
failure rates and test times of similar construction to the previous 10 systems, i.e., 
)t„(i) = KXj(i), where K is equal to 4.0 and X,t is also equal to 4.0. Each system also had 
its own random seeds for both DT and OT failures. Table 10 summarizes the four new 
systems. 


TABLE 10 

INITIAL DATA ON FOUR NEW "SIMILAR" SYSTEMS IN WHICH 
ACTUAL K = 4.0 AND ACTUAL POISSON FAILURE RATES X,t = 4.0 


System 

Actual DT 

DT test 

DT 

Actual OT 

OT test 

OT 

number 

failure 

time 

random 

failure 

time 

random 


rate 

(hours) 

seed 

rate 

(hours) 

seed 

11 

0.0003 

13,333.33 

12 

0.0012 

3,333.33 

112 

12 

0.0005 

8,000 

22 

0.002 

2,000 

122 

13 

0.0007 

5,714.29 

32 

0.0028 

1,428.57 1 

132 

14 

0.0009 

4,444.44 

42 

0.0036 

1,111.11 1 

142 
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One thousand repetitions of a simulated OT and DT were run on each of the four 
new systems, and then the exact same analysis as that done previously in section III.C.2 
was performed for each replication. Grand estimates appear in Tables 11 through 14 in 
the column entitled "Average value (n = 1000)" and were produced by averaging the 
estimates of the following from all 1000 simulated observations, or sub-estimates; 

1. The estimated OT failure rate derived from prediction, direct observation and 
weighted average; 

2. The estimated variance of the estimated OT failure rate derived from prediction, 
direct observation and weighted average; and 

3. The estimated standard error of the estimated OT failure rate derived from 
prediction, direct observation and weighted average. 

Then, in order to get an idea of the sample distribution of each of the above 
super-estimates, the sample variance, sample standard deviation and coefficient of 
variation® were computed for the 1000 sub-estimates for each system. Finally, the 
estimated mean-squared error, which is the sum of the square of the average of the biases 
with the average of the variances over all replications, i.e., 

MSB = (Average of the biases)^ +(Average of the variances), 
was calculated for each prediction method and placed in the last column of each 
table. Figure 3 describes the situation. 


®The coefficient of variation is simply the sample standard deviation divided by the sample average, and 
serves as a measure of the spread of the sample. 
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of estimating OT failure rate. 






1 . New System II. 


The statistics on new System 11 are summarized in Table 11 and a graphical 
summary is shown in Figure 4. Recall that actual (11) = 0.0012. 


TABLE 11 

NEW SYSTEM 11 - ESTIMATED OT FAILURE RATES AND OTHER 
RELEVANT STATISTICS (ACTUAL Xr,(l 1) = 0.0012) 



Sub- 

esdmates 

Average 

(n=1000) 

Sample 

Variance 

(n=1000) 

Sample Std 
Deviation 
(n=1000) 

Coeffof 

Variation 

Squared 

Prediction 


0.0012 

4.62x10-? 

0.00068 

57% 

4.46x10-? 

vrl„(ii)] 

4.46x10-'^ 

1.17x10-13 

3.43x10-? 

77% 

SE[i„(ll)] 

0.00062 

5.54x10-*! 

0.00024 1 

39% 

Direct 

Observation 

^Uii) 

0.0012 

3.73x10-? 

0.00061 

51% 

3.6x10-? 

va;(ii)] 

3.6x10-'^ 

3.36x10-1'* 

1.83x10-? 

51% 

SEfXUll)] 

0.00058 

2.87x10-? 

0.00017 

29% 

Weighted 

Average 

idi) 

0.0010 

2.23x10-? 

0.00047 

47% 

2.07x10-? 

vridi)] 

1.67x10-'? 

6.96x10-15 

8.34x10-*! 

50% 

SE[X,dl)] 

0.00039 

1.4x10-*! 

0.00012 

31% 


New Sysism 11 

■ Actual Failure Rale B Mean Estimated Failure Rate ■ Mean Estimated Std Error 



Prediction Direct Observation Weighted Average 


Figure 4 - Graphical summary of System 11 estimates. 
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b. New System 12. 


The statistics on new System 12 are summarized in Table 12 and a graphical 
summary is shown in Figure 5. Recall that actual (12) = 0.002. 


TABLE 12 

NEW SYSTEM 12 - ESTIMATED OT FAILURE RATES AND OTHER 
RELEVANT STATISTICS (ACTUAL X.r,(12) = 0.002) 



Sub- 

estimates 

Average 

(n=1000) 

Sample 

Variance 

(n=1000) 

Sample Std 
Deviation 
(n=1000) 

Coeffof 

Variation 

Squared 

Prediction 

X„(12) 

0.0020 

1.17x10-6 

0.0011 

55% 

1.23x10-6 

V[X.(12)] 

1.23x10-6 

7.51x10-13 

8.66x10-1 

70% 

SE[X,(12)] 

0.0010 

1.46x10-1 

0.00038 

38% 

Direct 

Observation 

^:(12) 

0.0020 

9.51x10-1 

0.00098 

49% 

9.99x10-1 

V[X;(12)] 

9.99x10-'^ 

2.38x10-13 

4.88x10-1 

49% 1 

SE[X:(12)] 

0.00096 

7.42x10-*! 

0.00027 

28% 

Weighted 

Average 

^.(12) 

0.0017 

5.91x10-1 

0.00077 

45% 

5.49x10-1 

V[i(12)] 

4.59x10-’^ 

S.BxlO-l'l 

2.26x10-1 

49% 

SE[1.(12)] 

0.00065 

3.71x10-*! 

0.00019 

29% 


NewSystsm 12 

■ Actual Failure Rate Q Mean Estimated Failure Rate ■ Mean Estimated Std Error 



Figure 5 - Graphical summary of System 12 estimates. 
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c. New System 13. 


The statistics on new System 13 are summarized in Table 13 and a graphical 
summary is shown in Figure 6. Recall that actual (13) = 0.0028. 


TABLE 13 

NEW SYSTEM 13 - ESTIMATED OT FAILURE RATES AND OTHER 
RELEVANT STATISTICS (ACTUAL A.r,(13) = 0.0028) 



Sub¬ 

estimates 

Average 

(n=1000) 

Sample 

Variance 

(n=1000) 

Sample Std 
Deviation 
(n=1000) 

Coeffof 

Variation 

Squared 

Error 

Prediction 


0.0028 

2.33x10-6 

0.0015 

54% 

2.41x10-6 

VU„(13)] 

2.41x10-6 

3.14x10-1^ 

1.77x10-6 

73% 

SE[X„(13)] 

0.0015 

2.96x10-7 

0.00054 

36% 

Direct 

Observation 

x:(i3) 

0.0028 

2.00x10-6 

0.0014 

50% 

1.96x10-6 

v[i:(i3)] 

1.96x10-6 

9.82x10-13 

9.91x10-7 

51% 

SE[X;(I3)] 

0.0013 

1.52x10-7 

0.00039 

30% 

Weighted 

XT13) 

0.0024 

1.24x10-6 

0.0011 

46% 


Average 

V[l„(13)] 

9.04x10-7 

2.11x10-13 

4.59x10-7 

51% 

1.06x10-6 


SE[i.(I3)] 

0.00091 

7.79x10-*^ 

0.00028 

31% 



■ Actual Failure Rale □ Mean Estimated Failure Rate ■ Mean Estimated Std Error 



Figure 6 - Graphical summary of System 13 estimates. 
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d. New System 14. 

The statistics on new System 14 are summarized in Table 14 and a graphical 
summary is shown in Figure 7. Recall that actual (14) = 0.0036. 


TABLE 14 

NEW SYSTEM 14 - ESTIMATED OT FAILURE RATES AND OTHER 
_ RELEVANT STATISTICS (ACTUAL A,r,(14) = 0.0036) 



Sub- 

estimaies 

Average 

(n=1000) 

Sample 

Variance 

(n=1000) 

Sample Std 
Deviation 
(n=1000) 

Coeffof 

Variation 

Squared 

Error 

Prediction 


0.0036 

4.07x10-6 

0.0020 

56% 

3.93x10-6 

V[X„(14)) 

3.93x10-6 

8.90x10-12 

2.98x10-6 

76% 

SE[X„(14)] 

0.0018 

5.18x10-2 

0.00072 

40% 

Direct 

Observation 


0.0036 

3.20x10-6 

0.0018 

50% 

3.21x10-6 

V[X:(14)] 

3.21x10-6 

2.59x10-12 

1.61x10-6 

50% 

SE[X:(14)] 

0.0017 

2.45x10-2 

0.00050 

29% 1 

Weighted 

i(14) 

0.0031 

2.05x10-6 

0.0014 

45% 


Average 

V(i(14)] 

1.47x10-6 

5.76x10-1^ 

7.59x10-2 

52% 

1.72x10-6 


SE[i.(14)] 

0.0012 

1.31x10-2 

0.00036 

30% 



New System 14 

■ Actual Failure Rate B Mean Estimated Failure Rate ■ Mean Estimated Std Error 



Figure 7 • Graphical summary of System 14 estimates. 
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D. RESULTS OF THE GAVER-JACOBS MODEL ANALYSIS 

Simulation results indicate that the predicted value method and the observed value 
method are unbiased and provide good estimates of the OT failure rate. The weighted 
average method, although it gives a more precise estimate because of its smaller mean 
squared error, is biased in that it tends to underestimate the actual failure rate. The 
observed value method gives slightly better precision than the predicted value method, 
since the mean squared error obtained with the observed value method tends to be slightly 
less than that associated with the predicted value method. Notice, however, that except 
for the very small price paid in terms of precision, the estimate obtained using the Gaver- 
Jacobs predicted value method seems to estimate the OT failure rate nearly as well as 
direct observation of the OT failure rate. This is significant since the Gaver-Jacobs 
estimate can be used prior to the start of operational testing (provided the constant-K 
model is nearly correct). 

The bottom line is that if the relationship between the DT failure rate and the OT 
failure rate proposed earlier does in fact exist, then the Gaver-Jacobs method of estimating 
the OT failure rate of a new system, A.„(i) = KA.j(i), given that system's DT failure rate as 
weU as DT and OT data on previous similar systems, seems to provide a reasonable 
estimate of the OT failure rate of the new system which performs nearly as well as direct 
observation of the OT failure rate through operational testing. Moreover, the Gaver- 
Jacobs estimate of the OT failure rate can be computed prior to the start of operational 
testing. 
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E. A HIERARCHICAL BAYES MODEL FOR ESTIMATION OF K 

It is idealistic, in the previous model, to assume that all systems have exactly the same 
value for K.® Rather, assume that each system i has its own K value, Kj, and that it is 
selected independently from a prior distribution. The gamma distribution, for example, 
turns out to be conjugate, so integrals can be explicitly computed. Its density function is: 

qK(K) = e-^'^^H^p. (3.22) 

From equations (3.5) and (3.6) it can be seen that the likelihood of Kj, up to irrelevant 
constants, is 


L,(K,;X/i);i-data) = e-''‘’‘-«‘-'‘'K;'’®[e-^'-«‘'-®(^d(i)r'‘’^°““^l- (3-23) 


To obtain the marginal posterior density of Kj, equation (3.22) is re-written, up to 
irrelevant constants, as follows; 


q^. (K,;i-data) = e"'' 


((p + X,(i)t„(i))K,r 
r(8-t-n„(i)) 


-(p + X,(i)t„(i)). (3.24) 


If the parameters p, 5 and ^j(i) were known, and the data available, i.e., n„(i) and 
t„(i), then one could compute the conditional expected value of Kj and use this as a point 
estimate. However there are other, more interesting ways to use the posterior. First, it is 

^The following explanation for the hierarchical Bayes model was developed by D. P. Gaver, Phofessor of 
Operations Research, Naval Postgraduate School, in handwritten notes to the author, July 1994. 
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jcessary to evaluate the various parameters p, 5 and A.j(i). To do this, the "empirical 


Bayes" method is used, meaning that the unknown K; term is integrated out of equation 
(3.23) to obtain 


L;= 


-(pK^r'p dK. fe-"' 




(3.25) 


The integral in the above expression can be reduced to 

^ J (3 26) 

Noting that the cumulative distribution function (CDF) of a gamma distribution with 
parameters (n„(i) + 5) and (A.j(i)t„(i) + p) has a known total area of 1.0, the following 
identity can be stated; 

^^'r(Va)"+S)°' ~J = 1. (3.27) 

Rearranging equation (3.26) to take advantage of (3.27), equation (3.26) can be reduced 


r(n„(i) + S)p^ 
r(8)(>.,(i)t„(i) + p)”-“>^ 


(3.28) 
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Therefore, equation (3.25) can be re-written as 




r(n„(i)-h5)p^ 


It is generally easier to work with the natural log of the Ukelihood function 


^. = ln(L,) = -A., (i)t<, (i) -i-(n, (i) + n„ (i))ln()L, (i)) -Hn(r(n„ (i) -t- 5)) 
-t- 5ln(p) -ln(r(8)) -(n„ (i) -t-8)ln(X.j (i)t„(i) -i- p) 


(3.30) 


The partial derivative of equation (3.30) with respect to (i) is 


(i) 


/■X 1 rx /-x —^d(i)to(i) 

n,(i)-A.j(i)tj(i) p _ 

A.,(i) ^ A.,(i)(A.,(i) + p) 


(3.31) 


Setting equation (3.31) equal to zero and solving for lj(i) gives the maximum likelihood 
estimator for (i). Unfortunately, however, the parameters S and p are unknown. 

One way to get at the parameters 5 and p is to obtain an initial estimate for (i) 
using the formula 


= (3.32) 

td(i) 

and then derive maximum likelihood estimators for 5 and p. One can then use these initial 
estimates for 5 and p in equation (3.31) to obtain an adjusted value for X,j(i) for each 
system. It may be possible to continue these iterations until 8, p, and (i) all stabilize. 
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Then, using 5 and p, as the parameters for the gamma prior distribution, it is possible t( 
estimate the mean Kj factor and its variance using the relationships 


Var[Ki] = ^. 
^ (P)^ 


Assuming data exists on more than one "similar" system, say 1 of them, then the log- 
likelihood function becomes, up to irrelevant constants 


^ = Sln(p)-i-ln -(n<,(i) + 5)ln(^d(i)to(i) + p) • (3-35) 


The ratio of gamma functions in equation (3.35) can be simplified by expanding the 
numerator and then canceling the denominator 

r(n„(i)-HS) (n„(i)-^S-l)(n„(i)-^S-2)-(n„(i)-^8-n„(i))^(8) 

r(8) ■ r(8) 

n„(i) = 0 (3.36) 

■|(n„(i) + 5-l)(n„(i) + 8-2)-(n„(i) + 6-n„(i)), n<,(i)>l. 


Then the natural log of equation (3.36) becomes 
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Rewriting equation (3.35) gives 


^ = g[6ln(p) + f(5)-(n„(i) + 8)ln(>L,(i)t„(i) + p)] (3.38) 

which can be further simplified to 

f = ^[5b(p)-(n<,(i)+8)ln(Xj(i)t„(i) + p)] + 5^f’(5). (3.39) 


It is now a simple matter to take the partial derivatives of equation (3.39) with respect to 
8 and p which, when set equal to zero, will enable one to solve for maximum likelihood 
estimators for the two parameters. Specifically 



(3.40) 


(3.41) 


(3.42) 


Finally, the preceding two equations, (3.40) and (3.42) along with (3.31) give 1 + 2 
equations with 1 + 2 unknowns which can be solved to obtain (i), 6 and p. 

So far, the author has not been able to solve the system of 1+2 equations and 1+2 
unknowns. Solution with I = 10 similar systems was attempted using the non-linear solver 
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included with Microsoft Excel 4.0 with solver options set as follows: Newton search 
method, using forward derivatives, with tangential estimates. Moreover, attempts to use 
the same methodology with only two "similar" systems; i.e., four equations and four 
unknowns, have also failed. The author believes that equations (3.31), (3.40) and (3.42), 
or at least some combination of the equations, are ill-conditioned in that they fail to 
converge at a common point in (1+2) space. An attempt was made to re-parameterize the 
equations, substituting the variables Xj(i), p and 8 with exponentials, however this 
attempt also failed to converge. 

While it appears the hierarchical Bayes methodology has promise, as well as a certain 
intuitive appeal, further research is called for in order to better condition the set of 
equations needed to use the model. Other issues to consider are whether or not the 
gamma prior is appropriate for this model and if so, how the set of equations needed to 
use the model can be simplified or re-parameterized so that reliable model parameter 
estimates can be obtained. 
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IV. SUMMARY, CONCLUSIONS AND RECOMMENDATIONS 


A. SUMMARY 

The puqjose of this thesis is to (1) investigate the behavior of a statistical procedure 
used to determine the length of time needed to test a continuously operating system when 
the system under test is assumed to have an exponential time to failure distribution, but 
when in fact the system does not have an exponential time to failure distribution and rather 
may have a time to failure distribution that more closely resembles a WeibuU distribution; 
(2) to examine the behavior of the Gaver-Jacobs estimator for the relationship between 
DT and OT failure rates using a Monte Carlo simulation; and (3) to introduce a 
hierarchical Bayes method for estimation of the relationship between the DT and OT 
failure rates. In the first part of the thesis, the parameters of the procedure are chosen 
based on the assumption of an exponential time to failure distribution. The behavior of the 
procedure is then investigated under the assumption that the time between failures are 
independent, having a WeibuU distribution with the same mean time to failure. The 
intention is to test whether or not use of the WeibuU distribution wUl result in fewer test 
hours to verify that a system under test has a mean time to faUure which meets or exceeds 
a minimum acceptable value at a given level of confidence, and also to test whether or not 
the OTAs are misstating the confidence they report in their operational testing results 
when the exponential distribution is used to model a system whose true failure rate 
function is that of a WeibuU distribution. In the second part of the thesis, the Gaver- 
Jacobs estimator is examined in terms of the consistency and variability of the OT failure 
rate derived using that estimator compared to other methods of establishing the OT faUure 
rate. Finally, in the third part of the thesis, a hierarchical Bayes model is proposed for the 
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ultimate purpose of being able to describe the mean and variance associated with a random 
relationship between DT and OT failure rates. 

B. CONCLUSIONS 

1. On comparison of the exponential and Weibull failure models. 

The first analysis suggests that operational testing agencies can, under certain 
circumstances, reduce the amount of time needed to demonstrate that the mean time to 
failure (MTTF) of a continuously operating system is at least as large as the minimum 
acceptable value (MAV) if the system's true underlying failure distribution is Weibull with 
shape parameter (3 > 1. However, the shape parameter is seldom, if ever, known or even 
estimable without expending valuable test resources, so it is unlikely that an OTA could 
actually save any time or resources with this methodology. Moreover, from the 
perspective of confidence levels, the analysis suggests that confidence in one's result is not 
sacrificed when the exponential distribution is used to model a system which has a non¬ 
constant failure rate, since it appears that use of the exponential distribution on a truly 
Weibull system seems to cause the OTA to overestimate the number of hours needed to 
test the MAV for the MTTF to a given confidence level, and also to underestimate the 
confidence it reports in its results under the incorrect assumption of an exponential 
distribution. However, for reasons mentioned previously, these results may not be true in 
general. Therefore, it is concluded that there is currently little to be gained by attempting 
to model the time between system failures with the Weibull distribution as previously 
suggested, particularly if very few observations of time to failure are available. 

2. On using the Gaver-Jacobs method to combine DT and OT failure data. 

The second analysis indicates that the Gaver-Jacobs method for estimating a 

system's OT failure rate given that system's DT failure rate as well as failure data on 
previous similar systems produces a reasonable estimate of the actual OT failure rate 
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(when the OT failure rate behaves as expected). In other words, if the relationship 
between the DT and OT failure rates for a certain system is such that it can in fact be 
estimated from previous similar systems, then the Gaver-Jacobs method gives a good 
estimate of the OT failure rate given that system's DT failure rate. In fact, the estimate 
provided by the Gaver-Jacobs method appears to be nearly as good as direct observation 
of the OT failure rate obtained through actual testing. Thus, if it can be shown that 
inferences about the failure rate of a system can in fact be made on the basis of previous 
similar systems, then the Gaver-Jacobs method of predicting a new system's OT failure 
rate from previous DT failure data has some potential utility. In particular, this OT failure 
rate prediction could be used in ways which would aid decisions as to when to begin OT 
as well as to augment follow-on OT in an effort to reduce needed testing resources. 

3. On the hierarchical Bayes Model for estimating K. 

The third analysis is incomplete. It was hoped that the hierarchical Bayes method 
given would provide a means by which to estimate K, the relationship between the DT and 
OT failure rates, when K is assumed to be a random variable with some prior distribution. 
Regretfully, the equations derived for the maximum likelihood estimators for the 
parameters of the gamma prior have proved to be difficult to solve satisfactorily. The 
author was unable to find a meaningful solution to the set of 1+2 equations and 1+2 
unknowns, where 1 represents the number of "similar" systems for which previous DT and 
OT data exist. If such a solution were available it could be used to predict future failure 
rates by means of a Bayesian calculation. Investigation of such a procedure is left for a 
later time. 

C. RECOMMENDATIONS FOR FURTHER STUDY 

Concerning the notion of modeling the failure distribution of a continuous system 
using the WeibuU distribution, OTAs could clearly benefit from being able to use this 
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approach when appropriate. If the time between failures for a system has a WeibuU 
distribution, then proper choice of that system's WeibuU shape parameter would probably 
enable the OTA to reduce the required test resources while maintaining a specified level of 
confidence in the result of their test. Unfortunately, there is very httle data in existence 
which can assist analysts in estimating the shape parameter of a system unless it is assumed 
that shape parameter information can be "borrowed" from similar systems. Moreover, it is 
seldom worth the effort needed to submit the system to special testing just to determine 
it’s shape parameter. FinaUy, there is some doubt that the necessary relationship between 
expected time to failure and the numbers of failures observed in a fixed time exists 
between the exponential and WeibuU distributions. Therefore, one area for further 
research would be in devising methods by which analysts could estimate the shape 
parameter of a system prior to testing; see the "borrowing information" idea suggested 
above. Another area of further research, which would be useful in the area previously 
mentioned, would be in creating a database of system types and the associated WeibuU 
shape parameters of their inter-failure distribution. A third area of investigation would be 
in identifying those combinations of exponential and WeibuU parameters in which a 
suitable relationship exists between the expected time to failure and the numbers of 
failures observed in a fixed time. 

Concerning the combining of DT and OT failure data, the abiUty to make inferences 
about the OT faUure rate of a system prior to testing would certainly benefit the OTAs by 
enabling them to defer testing until they are reasonably sure the system wiU pass the test. 
Moreover, the ability to augment future OT failure data with pseudo-data obtained from 
prior DT failure data would enable the OTAs to reduce the amount of testing resources 
needed to verify reliabUity thresholds. To date, the notion of inferring the relationship 
between DT failure data and OT failure data based on previous similar systems is untested. 
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Therefore, further areas of research in this field include the investigation of actual data on 
several similar systems to determine the feasibility of detecting any relationship that might 
exist between them and also to investigate the nature of the relationship which is found. 

Finally, concerning the estimation of the K factor using a hierarchical Bayes 
methodology, further research is needed in order to develop a more stable set of equations 
for the parameters of the gamma prior. Moreover, one could consider whether or not a 
gamma prior is appropriate for this model, or whether some other prior distribution might 
be more appropriate. 


50 



APPENDIX A. ALTERNATIVE METHOD OF DERIVING TIME TO TEST 

There is an alternative way to determine, for a given level of confidence, the length of 
time, t, needed to test a continuous-type system in order to demonstrate that the true 
MTTF of the system is at least as large as the MAV.if’ 

Let the MAV be denoted by 0, the 100a percent lower confidence limit for 7 , the true 
MTTF. The usual test plan is to test for a length of time, t, and observe F, the random 
number of failures to occur in (0, t]. The planned test time, t, is chosen so that if one 
failure is observed, the computed lower confidence limit, 7 ^. will be equal to the MAV, 0. 

The lower confidence limit, is the solution for 7 in the equation 


]-a = P(T>t) = P(N(t)<F) = ;£- 




where T is the time to the (F+l)st failure. 

The following relationship can be used to solve for 7 in equation (A.l) 




where xl^ denotes a chi-square random variable with 2a degrees of freedom. Equation 
(A.l) then becomes 


author is deeply indebted to W. Max Woods, Professor of Operations Research, Naval Postgraduate 
School, for his guidance in the derivation of this alternative method. 
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(A.3) 


1 a - P % 2 (f+i) - 


which implies that 


2t 2 2t 

- = Y=-5- (A-4) 

I A.a.2(F+l) 

where X^ jcf+i) o^th percentile of the chi-square distribution with 2 (F-f 1) degrees of 
freedom. Since the solution for 7 is and the MAV is 0, 

9 = Yl=^^- (A-5) 

Xa2(F+l) 


Therefore 


(A.6) 


If F = 1 and a = 0.8, then = 5.999 and thus, from equation (A.6), t = 30, as 
shown earlier in section n.A. Note that equation (A.6) holds for aU positive values of 0 
and F, making this a very simple way to solve for t given any values of 0 and F. 
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APPENDIX B. CODE FOR RANDOM POISSON ARRIVAL GENERATOR 


In this appendix, the Poisson failure generator discussed in Chapter IH is given. The 
program was used to simulate the number of failures during a given phase of testing for a 
given exponential failure rate and a given test duration. The program was written in 
Pascal and executed using Turbo Pascal 6.0 (Borland, 1990), a commercial Pascal 
package. 

Program Generate_Failures; 
uses TOOLBOX, CRT; 

{ Author: Timothy P. Anderson, Naval Postgraduate School 
Date: 27 March 1994 

Purpose: To generate numbers of failures that occur during a pre specified 
time, TIME, when the underlying failure distribution is 
exponential with rate LAMBDA. This program is written in partial 
support of a thesis by the author. 


var 

Lambda, Time, X, CumTime, Rand : real; 

I, N, Count: integer; 

View: boolean; 

OutFile: string; 

Ch : char; 

begin 

CRT.Clrscr; 

writeln ('Welcome to Program Generate_Failures... Follow instructions...'); 
writeln; writeln; 

write ('Specify the random seed:> '); readln (randseed); writeln; 
write ('Specify the exponential failure rate, Lambda:> '); 

Lambda := TOOLBOX.Get_Real; writeln; 
write ('Specify the testing time allowed:> '); 

Time:= TOOLBOX.Get_Real; writeln; 
write ('Specify number of iterations to run:> '); 

Count: = TOOLBOX.GetJnteger; writeln; 

writeln ('You may view output on screen or send it to a text file.'); 
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write ('Do you want to view it on the screen? (Y/N):> '); 

View := TOOLBOX.Get_Answer; writeln; writeln; writeln; 
if View = False then begin 
write ('Specify the output file (including path):> '); 
read (OutFile); readln; 

SYSTEM.Assign (Output, OutFile); 
rewrite (Output); 
end; {if} 

for I := 1 to Count do begin 
X:=0; 

N :=0; 

CumTime := 0; 
repeat 

Rand ;= Random; 

X := -(1/Lambda)*ln(l - Rand); 
if (X + CumTime) > Time then begin 
writeln (N); 

CumTime := CumTime + X; 
end 

else begin 
N :=N + 1; 

CumTime ;= CumTime + X; 
end; {if} 

until CumTime > Time; 
end; {for} 

SYSTEM.Close (Output); 

SYSTEM. Assign (Output, 'con'); 
rewrite (Output); 

if View = False then begin 

writeln ('Output sent to ', OutFile,'... Press any key to exit.'); 
Ch := Readkey; 
end 

else begin 

writeln ('Press any key to exit.'); 

Ch := Readkey; 
end; 

end. 
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